home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / pcl / src-16f.lha / code / irrat.lisp < prev    next >
Encoding:
Text File  |  1992-05-30  |  12.5 KB  |  395 lines
  1. ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
  2. ;;;
  3. ;;; **********************************************************************
  4. ;;; This code was written as part of the CMU Common Lisp project at
  5. ;;; Carnegie Mellon University, and has been placed in the public domain.
  6. ;;; If you want to use this code or any part of CMU Common Lisp, please contact
  7. ;;; Scott Fahlman or slisp-group@cs.cmu.edu.
  8. ;;;
  9. (ext:file-comment
  10.   "$Header: irrat.lisp,v 1.10 92/02/14 23:45:06 wlott Locked $")
  11. ;;;
  12. ;;; **********************************************************************
  13. ;;;
  14. ;;; This file contains all the irrational functions.  Actually, most of the
  15. ;;; work is done by calling out to C...
  16. ;;;
  17. ;;; Author: William Lott.
  18. ;;; 
  19.  
  20. (in-package "KERNEL")
  21.  
  22.  
  23. ;;;; Random constants, utility functions, and macros.
  24.  
  25. (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
  26. ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
  27.  
  28. ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
  29. ;;; call, and saves number consing to boot.
  30. ;;;
  31. (defmacro def-math-rtn (name num-args)
  32.   (let ((function (intern (concatenate 'simple-string
  33.                        "%"
  34.                        (string-upcase name)))))
  35.     `(progn
  36.        (proclaim '(inline ,function))
  37.        (export ',function)
  38.        (alien:def-alien-routine (,name ,function) double-float
  39.      ,@(let ((results nil))
  40.          (dotimes (i num-args (nreverse results))
  41.            (push (list (intern (format nil "ARG-~D" i))
  42.                'double-float)
  43.              results)))))))
  44.  
  45. (eval-when (compile load eval)
  46.  
  47. (defun handle-reals (function var)
  48.   `((((foreach fixnum single-float bignum ratio))
  49.      (coerce (,function (coerce ,var 'double-float)) 'single-float))
  50.     ((double-float)
  51.      (,function ,var))))
  52.  
  53. ); eval-when (compile load eval)
  54.  
  55.  
  56. ;;;; Stubs for the Unix math library.
  57.  
  58. ;;; Please refer to the Unix man pages for details about these routines.
  59.  
  60. ;;; Trigonometric.
  61. (def-math-rtn "sin" 1)
  62. (def-math-rtn "cos" 1)
  63. (def-math-rtn "tan" 1)
  64. (def-math-rtn "asin" 1)
  65. (def-math-rtn "acos" 1)
  66. (def-math-rtn "atan" 1)
  67. (def-math-rtn "atan2" 2)
  68. (def-math-rtn "sinh" 1)
  69. (def-math-rtn "cosh" 1)
  70. (def-math-rtn "tanh" 1)
  71. (def-math-rtn "asinh" 1)
  72. (def-math-rtn "acosh" 1)
  73. (def-math-rtn "atanh" 1)
  74.  
  75. ;;; Exponential and Logarithmic.
  76. (def-math-rtn "exp" 1)
  77. (def-math-rtn "expm1" 1)
  78. (def-math-rtn "log" 1)
  79. (def-math-rtn "log10" 1)
  80. (def-math-rtn "log1p" 1)
  81. (def-math-rtn "pow" 2)
  82. (def-math-rtn "cbrt" 1)
  83. (def-math-rtn "sqrt" 1)
  84. (def-math-rtn "hypot" 2)
  85.  
  86.  
  87. ;;;; Power functions.
  88.  
  89. (defun exp (number)
  90.   "Return e raised to the power NUMBER."
  91.   (number-dispatch ((number number))
  92.     (handle-reals %exp number)
  93.     ((complex)
  94.      (* (exp (realpart number))
  95.     (cis (imagpart number))))))
  96.  
  97. ;;; INTEXP -- Handle the rational base, integer power case.
  98.  
  99. (defparameter *intexp-maximum-exponent* 10000)
  100.  
  101. ;;; This function precisely calculates base raised to an integral power.  It
  102. ;;; separates the cases by the sign of power, for efficiency reasons, as powers
  103. ;;; can be calculated more efficiently if power is a positive integer.  Values
  104. ;;; of power are calculated as positive integers, and inverted if negative.
  105. ;;;
  106. (defun intexp (base power)
  107.   (when (> (abs power) *intexp-maximum-exponent*)
  108.     (cerror "Continue with calculation."
  109.         "The absolute value of ~S exceeds ~S."
  110.         power '*intexp-maximum-exponent* base power))
  111.   (cond ((minusp power)
  112.      (/ (intexp base (- power))))
  113.     ((eql base 2)
  114.      (ash 1 power))
  115.     (t
  116.      (do ((nextn (ash power -1) (ash power -1))
  117.           (total (if (oddp power) base 1)
  118.              (if (oddp power) (* base total) total)))
  119.          ((zerop nextn) total)
  120.        (setq base (* base base))
  121.        (setq power nextn)))))
  122.  
  123.  
  124. ;;; EXPT  --  Public
  125. ;;;
  126. ;;;    If an integer power of a rational, use INTEXP above.  Otherwise, do
  127. ;;; floating point stuff.  If both args are real, we try %POW right off,
  128. ;;; assuming it will return 0 if the result may be complex.  If so, we call
  129. ;;; COMPLEX-POW which directly computes the complex result.  We also separate
  130. ;;; the complex-real and real-complex cases from the general complex case.
  131. ;;;
  132. (defun expt (base power)
  133.   "Returns BASE raised to the POWER."
  134.   (if (zerop power)
  135.       (1+ (* base power))
  136.       (labels ((real-expt (base power rtype)
  137.          (let* ((fbase (coerce base 'double-float))
  138.             (fpower (coerce power 'double-float))
  139.             (res (coerce (%pow fbase fpower) rtype)))
  140.            (if (and (zerop res) (minusp fbase))
  141.                (multiple-value-bind (re im)
  142.                         (complex-pow fbase fpower)
  143.              (%make-complex (coerce re rtype) (coerce im rtype)))
  144.                res)))
  145.            (complex-pow (fbase fpower)
  146.          (let ((pow (%pow (- fbase) fpower))
  147.                (fpower*pi (* fpower pi)))
  148.            (values (* pow (%cos fpower*pi))
  149.                (* pow (%sin fpower*pi))))))
  150.     (declare (inline real-expt))
  151.     (number-dispatch ((base number) (power number))
  152.       (((foreach fixnum (or bignum ratio) (complex rational)) integer)
  153.        (intexp base power))
  154.       (((foreach single-float double-float) rational)
  155.        (real-expt base power '(dispatch-type base)))
  156.       (((foreach fixnum (or bignum ratio) single-float)
  157.         (foreach ratio single-float))
  158.        (real-expt base power 'single-float))
  159.       (((foreach fixnum (or bignum ratio) single-float double-float)
  160.         double-float)
  161.        (real-expt base power 'double-float))
  162.       ((double-float single-float)
  163.        (real-expt base power 'double-float))
  164.       (((foreach (complex rational) (complex float)) rational)
  165.        (* (expt (abs base) power)
  166.           (cis (* power (phase base)))))
  167.       (((foreach fixnum (or bignum ratio) single-float double-float)
  168.         complex)
  169.        (if (minusp base)
  170.            (/ (exp (* power (truly-the float (log (- base))))))
  171.            (exp (* power (truly-the float (log base))))))
  172.       (((foreach (complex float) (complex rational)) complex)
  173.        (exp (* power (log base))))))))
  174.  
  175. (defun log (number &optional (base nil base-p))
  176.   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
  177.   (if base-p
  178.       (/ (log number) (log base))
  179.       (number-dispatch ((number number))
  180.     (((foreach fixnum bignum ratio single-float))
  181.      (if (minusp number)
  182.          (complex (log (- number)) (coerce pi 'single-float))
  183.          (coerce (%log (coerce number 'double-float)) 'single-float)))
  184.     ((double-float)
  185.      (if (minusp number)
  186.          (complex (log (- number)) (coerce pi 'double-float))
  187.          (%log number)))
  188.     ((complex) (complex (log (abs number)) (phase number))))))
  189.  
  190. (defun sqrt (number)
  191.   "Return the square root of NUMBER."
  192.   (number-dispatch ((number number))
  193.     (((foreach fixnum bignum ratio single-float))
  194.      (if (minusp number)
  195.      (exp (/ (log number) 2))
  196.      (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
  197.     ((double-float)
  198.      (if (minusp number)
  199.      (exp (/ (log number) 2))
  200.      (%sqrt number)))
  201.     ((complex) (exp (/ (log number) 2)))))
  202.  
  203. ;;; ISQRT:  Integer square root - isqrt(n)**2 <= n
  204. ;;; Upper and lower bounds on the result are estimated using integer-length.
  205. ;;; On each iteration, one of the bounds is replaced by their mean.
  206. ;;; The lower bound is returned when the bounds meet or differ by only 1.
  207. ;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.
  208.  
  209. (defun isqrt (n)
  210.   "Returns the root of the nearest integer less than
  211.    n which is a perfect square."
  212.   (if (and (integerp n) (not (minusp n)))
  213.       (do* ((lg (integer-length n))
  214.         (lo (ash 1 (ash (1- lg) -1)))
  215.         (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.
  216.        ((<= (1- hi) lo) lo)
  217.     (let ((mid (ash (+ lo hi) -1)))
  218.       (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))
  219.       (error "Isqrt: ~S argument must be a nonnegative integer" n)))
  220.  
  221.  
  222.  
  223. ;;;; Trigonometic and Related Functions
  224.  
  225. (defun abs (number)
  226.   "Returns the absolute value of the number."
  227.   (number-dispatch ((number number))
  228.     (((foreach single-float double-float fixnum rational))
  229.      (abs number))
  230.     ((complex)
  231.      (let ((rx (realpart number))
  232.        (ix (imagpart number)))
  233.        (etypecase rx
  234.      (rational
  235.       (sqrt (+ (* rx rx) (* ix ix))))
  236.      (single-float
  237.       (coerce (%hypot (coerce rx 'double-float)
  238.               (coerce ix 'double-float))
  239.           'single-float))
  240.      (double-float
  241.       (%hypot rx ix)))))))
  242.  
  243. (defun phase (number)
  244.   "Returns the angle part of the polar representation of a complex number.
  245.   For complex numbers, this is (atan (imagpart number) (realpart number)).
  246.   For non-complex positive numbers, this is 0.  For non-complex negative
  247.   numbers this is PI."
  248.   (etypecase number
  249.     ((or rational single-float)
  250.      (if (minusp number)
  251.      (coerce pi 'single-float)
  252.      0.0f0))
  253.     (double-float
  254.      (if (minusp number)
  255.      (coerce pi 'double-float)
  256.      0.0d0))
  257.     (complex
  258.      (atan (imagpart number) (realpart number)))))
  259.  
  260.  
  261. (defun sin (number)  
  262.   "Return the sine of NUMBER."
  263.   (number-dispatch ((number number))
  264.     (handle-reals %sin number)
  265.     ((complex)
  266.      (let ((x (realpart number))
  267.        (y (imagpart number)))
  268.        (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
  269.  
  270. (defun cos (number)
  271.   "Return the cosine of NUMBER."
  272.   (number-dispatch ((number number))
  273.     (handle-reals %cos number)
  274.     ((complex)
  275.      (let ((x (realpart number))
  276.        (y (imagpart number)))
  277.        (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
  278.  
  279. (defun tan (number)
  280.   "Return the tangent of NUMBER."
  281.   (number-dispatch ((number number))
  282.     (handle-reals %tan number)
  283.     ((complex)
  284.      (let* ((num (sin number))
  285.         (denom (cos number)))
  286.        (if (zerop denom) (error "~S undefined tangent." number)
  287.        (/ num denom))))))
  288.  
  289. (defun cis (theta)
  290.   "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
  291.   (if (complexp theta)
  292.       (error "Argument to CIS is complex: ~S" theta)
  293.       (complex (cos theta) (sin theta))))
  294.  
  295. (proclaim '(inline mult-by-i))
  296. (defun mult-by-i (number)
  297.   (complex (imagpart number)
  298.        (- (realpart number))))
  299.  
  300. (defun complex-asin (number)
  301.   (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
  302.  
  303. (defun asin (number)
  304.   "Return the arc sine of NUMBER."
  305.   (number-dispatch ((number number))
  306.     ((rational)
  307.      (if (or (> number 1) (< number -1))
  308.      (complex-asin number)
  309.      (coerce (%asin (coerce number 'double-float)) 'single-float)))
  310.     (((foreach single-float double-float))
  311.      (if (or (> number (coerce 1 '(dispatch-type number)))
  312.          (< number (coerce -1 '(dispatch-type number))))
  313.      (complex-asin number)
  314.      (coerce (%asin (coerce number 'double-float))
  315.          '(dispatch-type number))))
  316.     ((complex)
  317.      (complex-asin number))))
  318.  
  319. (defun complex-acos (number)
  320.   (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
  321.  
  322. (defun acos (number)
  323.   "Return the arc cosine of NUMBER."
  324.   (number-dispatch ((number number))
  325.     ((rational)
  326.      (if (or (> number 1) (< number -1))
  327.      (complex-acos number)
  328.      (coerce (%acos (coerce number 'double-float)) 'single-float)))
  329.     (((foreach single-float double-float))
  330.      (if (or (> number (coerce 1 '(dispatch-type number)))
  331.          (< number (coerce -1 '(dispatch-type number))))
  332.      (complex-acos number)
  333.      (coerce (%acos (coerce number 'double-float))
  334.          '(dispatch-type number))))
  335.     ((complex)
  336.      (complex-acos number))))
  337.  
  338.  
  339. (defun atan (y &optional (x nil xp))
  340.   "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
  341.   (if xp
  342.       (if (and (zerop x) (zerop y))
  343.       (if (plusp (float-sign (float x)))
  344.           y
  345.           (if (minusp (float-sign (float y)))
  346.           (- pi)
  347.           pi))
  348.       (number-dispatch ((y real) (x real))
  349.         (((foreach fixnum bignum ratio single-float)
  350.           (foreach fixnum bignum ratio single-float))
  351.          (coerce (%atan2 (coerce y 'double-float)
  352.                  (coerce x 'double-float))
  353.              'single-float))
  354.         ((double-float (foreach fixnum bignum ratio single-float))
  355.          (%atan2 y (coerce x 'double-float)))
  356.         (((foreach fixnum bignum ratio single-float double-float)
  357.           double-float)
  358.          (%atan2 (coerce y 'double-float) x))))
  359.       (number-dispatch ((y number))
  360.     (handle-reals %atan y)
  361.     ((complex)
  362.      (let ((im (imagpart y))
  363.            (re (realpart y)))
  364.        (/ (- (log (complex (- 1 im) re))
  365.          (log (complex (+ 1 im) (- re))))
  366.           (complex 0 2)))))))
  367.  
  368.  
  369. (defun sinh (number)
  370.   "Return the hyperbolic sine of NUMBER."
  371.   (/ (- (exp number) (exp (- number))) 2))
  372.  
  373. (defun cosh (number)
  374.   "Return the hyperbolic cosine of NUMBER."
  375.   (/ (+ (exp number) (exp (- number))) 2))
  376.  
  377. (defun tanh (number)
  378.   "Return the hyperbolic tangent of NUMBER."
  379.   (/ (- (exp number) (exp (- number)))
  380.      (+ (exp number) (exp (- number)))))
  381.  
  382.  
  383. (defun asinh (number)
  384.   "Return the hyperbolic arc sine of NUMBER."
  385.   (log (+ number (sqrt (1+ (* number number))))))
  386.  
  387. (defun acosh (number)
  388.   "Return the hyperbolic arc cosine of NUMBER."
  389.   (log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
  390.  
  391. (defun atanh (number)
  392.   "Return the hyperbolic arc tangent of NUMBER."
  393.   (log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
  394.  
  395.